packages <- c('maptools', 'sf', 'raster', 'spatstat', 'tmap', 'tidyverse','plotly', 'ggthemes')
for(p in packages){
if(!require(p, character.only=T)){
install.packages(p)
}
library(p,character.only = T)
}
CHAS <- read_rds("data/rds/CHAS.rds")
childcare <- read_rds("data/rds/childcare.rds")
mpsz_sf <- st_read(dsn = "data/shapefile", layer = "MP14_SUBZONE_WEB_PL")
Reading layer `MP14_SUBZONE_WEB_PL' from data source
`C:\yiling-yu\IS415_Blog\_posts\2021-09-13-in-class-exercise-5\data\shapefile'
using driver `ESRI Shapefile'
Simple feature collection with 323 features and 15 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 2667.538 ymin: 15748.72 xmax: 56396.44 ymax: 50256.33
Projected CRS: SVY21
glimpse(CHAS)
Rows: 1,192
Columns: 17
$ NAME <chr> "FONG CLINIC", "FRASER MEDICAL CENTR~
$ DESCRIPTION <chr> "CDMP,CHAS", "CDMP,CHAS", "CDMP,CHAS~
$ ADDRESSBLOCKHOUSENUMBER <chr> "162", "78A", "1", "982", "163", "30~
$ ADDRESSFLOORNUMBER <chr> "01", "01", "01", "01", "01", "01", ~
$ ADDRESSPOSTALCODE <chr> "140162", "101078", "150001", "53098~
$ ADDRESSSTREETNAME <chr> "MEI LING ST", "TELOK BLANGAH ST 32"~
$ ADDRESSUNITNUMBER <chr> "361", "07", "4524", "03", "426", "8~
$ X_COORDINATE <chr> "24609.88367388", "25288.1926698", "~
$ Y_COORDINATE <chr> "30516.49536544", "28423.68094962", ~
$ HCI_CODE <chr> "9403524", "16C0164", "9400550", "94~
$ LICENCE_TYPE <chr> "MC", "MC", "MC", "MC", "MC", "MC", ~
$ HCI_TEL <chr> "64745993", "62733603", "62726628", ~
$ ADDR_TYPE <chr> "A", "A", "A", "A", "A", "A", "A", "~
$ Lat <chr> "1.2922546717017", "1.2733280593593"~
$ Lng <chr> "103.802856952292", "103.80895202390~
$ ICON_NAME <chr> "OM2_MOHicon_CHASclinic_20px.png", "~
$ ADDRESSBUILDINGNAME <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, ~
glimpse(childcare)
Rows: 1,545
Columns: 7
$ NAME <chr> "AVERBEL CHILD DEVELOPMENT CENTRE PTE LTD"~
$ DESCRIPTION <chr> "Child Care Services", "Child Care Service~
$ ADDRESSPOSTALCODE <chr> "760742", "159053", "556912", "569139", "4~
$ ADDRESSSTREETNAME <chr> "742, YISHUN AVENUE 5, #01 - 470, SINGAPOR~
$ Lat <chr> "1.42972006815634", "1.28668026511187", "1~
$ Lng <chr> "103.833109448847", "103.813766346615", "1~
$ ICON_NAME <chr> "onemap-fc-childcare.png", "onemap-fc-chil~
mpsz_sf <- st_set_crs(mpsz_sf, 3414)
st_crs(mpsz_sf)
Coordinate Reference System:
User input: EPSG:3414
wkt:
PROJCRS["SVY21 / Singapore TM",
BASEGEOGCRS["SVY21",
DATUM["SVY21",
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4757]],
CONVERSION["Singapore Transverse Mercator",
METHOD["Transverse Mercator",
ID["EPSG",9807]],
PARAMETER["Latitude of natural origin",1.36666666666667,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",103.833333333333,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["Scale factor at natural origin",1,
SCALEUNIT["unity",1],
ID["EPSG",8805]],
PARAMETER["False easting",28001.642,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["False northing",38744.572,
LENGTHUNIT["metre",1],
ID["EPSG",8807]]],
CS[Cartesian,2],
AXIS["northing (N)",north,
ORDER[1],
LENGTHUNIT["metre",1]],
AXIS["easting (E)",east,
ORDER[2],
LENGTHUNIT["metre",1]],
USAGE[
SCOPE["Cadastre, engineering survey, topographic mapping."],
AREA["Singapore - onshore and offshore."],
BBOX[1.13,103.59,1.47,104.07]],
ID["EPSG",3414]]
CHAS_sf <- st_as_sf(CHAS,
coords = c("X_COORDINATE",
"Y_COORDINATE"),
crs=3414)
childcare_sf <- st_as_sf(childcare,
coords = c("Lng",
"Lat"),
crs=4326) %>%
st_transform(crs = 3414)
#plot the graph to have a look
tmap_mode("view")
tm_shape(childcare_sf) +
tm_dots(alpha=0.4,
col= "blue",
size = 0.05) +
tm_shape(CHAS_sf) +
tm_dots(alpha=0.4,
col= "red",
size = 0.05)
childcare <- as_Spatial(childcare_sf)
CHAS <- as_Spatial(CHAS_sf)
mpsz <- as_Spatial(mpsz_sf)
childcare_sp <- as(childcare, "SpatialPoints")
CHAS_sp <- as(CHAS, "SpatialPoints")
mpsz_sp <- as(mpsz, "SpatialPolygons")
childcare_ppp <- as(childcare_sp,"ppp")
CHAS_ppp <- as(CHAS_sp,"ppp")
childcare_ppp_jit <- rjitter(childcare_ppp,
retry = TRUE,
nsim = 1,
drop = TRUE)
any(duplicated(childcare_ppp_jit))
[1] FALSE
CHAS_ppp_jit <- rjitter(CHAS_ppp,
retry = TRUE,
nsim = 1,
drop = TRUE)
any(duplicated(CHAS_ppp_jit))
[1] FALSE
pg <- mpsz[mpsz@data$PLN_AREA_N=="PUNGGOL",]
pg_sp <- as(pg, "SpatialPolygons")
pg_owin <- as(pg_sp, "owin")
childcare_pg <- childcare_ppp_jit[pg_owin]
CHAS_pg <- CHAS_ppp_jit[pg_owin]
plot(childcare_ppp_jit)
plot(childcare_pg)

L_childcare <- envelope(childcare_pg,
Lest,
nsim = 99,
rank = 1,
gloval = TRUE)
Generating 99 simulations of CSR ...
1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35,
36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70,
71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99.
Done.
L_CHAS <- envelope(CHAS_pg,
Lest,
nsim = 99,
rank = 1,
gloval = TRUE)
Generating 99 simulations of CSR ...
1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35,
36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70,
71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99.
Done.
title <- "Pairwise Distance: L function"
Lcsr_df <- as.data.frame(L_childcare)
colour=c("#0D657D","#ee770d","#D3D3D3")
csr_plot <- ggplot(Lcsr_df, aes(r, obs-r))+
# plot observed value
geom_line(colour=c("#4d4d4d"))+
geom_line(aes(r,theo-r), colour="red", linetype = "dashed")+
# plot simulation envelopes
geom_ribbon(aes(ymin=lo-r,ymax=hi-r),alpha=0.1, colour=c("#91bfdb")) +
xlab("Distance r (m)") +
ylab("L(r)-r") +
geom_rug(data=Lcsr_df[Lcsr_df$obs > Lcsr_df$hi,], sides="b", colour=colour[1]) +
geom_rug(data=Lcsr_df[Lcsr_df$obs < Lcsr_df$lo,], sides="b", colour=colour[2]) +
geom_rug(data=Lcsr_df[Lcsr_df$obs >= Lcsr_df$lo & Lcsr_df$obs <= Lcsr_df$hi,], sides="b", color=colour[3]) +
theme_tufte()+
ggtitle(title)
text1<-"Significant clustering"
text2<-"Significant segregation"
text3<-"Not significant clustering/segregation"
# the below conditional statement is required to ensure that the labels (text1/2/3) are assigned to the correct traces
if (nrow(Lcsr_df[Lcsr_df$obs > Lcsr_df$hi,])==0){
if (nrow(Lcsr_df[Lcsr_df$obs < Lcsr_df$lo,])==0){
ggplotly(csr_plot, dynamicTicks=T) %>%
style(text = text3, traces = 4) %>%
rangeslider()
}else if (nrow(Lcsr_df[Lcsr_df$obs >= Lcsr_df$lo & Lcsr_df$obs <= Lcsr_df$hi,])==0){
ggplotly(csr_plot, dynamicTicks=T) %>%
style(text = text2, traces = 4) %>%
rangeslider()
}else {
ggplotly(csr_plot, dynamicTicks=T) %>%
style(text = text2, traces = 4) %>%
style(text = text3, traces = 5) %>%
rangeslider()
}
} else if (nrow(Lcsr_df[Lcsr_df$obs < Lcsr_df$lo,])==0){
if (nrow(Lcsr_df[Lcsr_df$obs >= Lcsr_df$lo & Lcsr_df$obs <= Lcsr_df$hi,])==0){
ggplotly(csr_plot, dynamicTicks=T) %>%
style(text = text1, traces = 4) %>%
rangeslider()
} else{
ggplotly(csr_plot, dynamicTicks=T) %>%
style(text = text1, traces = 4) %>%
style(text = text3, traces = 5) %>%
rangeslider()
}
} else{
ggplotly(csr_plot, dynamicTicks=T) %>%
style(text = text1, traces = 4) %>%
style(text = text2, traces = 5) %>%
style(text = text3, traces = 6) %>%
rangeslider()
}
title <- "Pairwise Distance: L function"
Lcsr_df <- as.data.frame(L_CHAS)
colour=c("#0D657D","#ee770d","#D3D3D3")
csr_plot <- ggplot(Lcsr_df, aes(r, obs-r))+
# plot observed value
geom_line(colour=c("#4d4d4d"))+
geom_line(aes(r,theo-r), colour="red", linetype = "dashed")+
# plot simulation envelopes
geom_ribbon(aes(ymin=lo-r,ymax=hi-r),alpha=0.1, colour=c("#91bfdb")) +
xlab("Distance r (m)") +
ylab("L(r)-r") +
geom_rug(data=Lcsr_df[Lcsr_df$obs > Lcsr_df$hi,], sides="b", colour=colour[1]) +
geom_rug(data=Lcsr_df[Lcsr_df$obs < Lcsr_df$lo,], sides="b", colour=colour[2]) +
geom_rug(data=Lcsr_df[Lcsr_df$obs >= Lcsr_df$lo & Lcsr_df$obs <= Lcsr_df$hi,], sides="b", color=colour[3]) +
theme_tufte()+
ggtitle(title)
text1<-"Significant clustering"
text2<-"Significant segregation"
text3<-"Not significant clustering/segregation"
# the below conditional statement is required to ensure that the labels (text1/2/3) are assigned to the correct traces
if (nrow(Lcsr_df[Lcsr_df$obs > Lcsr_df$hi,])==0){
if (nrow(Lcsr_df[Lcsr_df$obs < Lcsr_df$lo,])==0){
ggplotly(csr_plot, dynamicTicks=T) %>%
style(text = text3, traces = 4) %>%
rangeslider()
}else if (nrow(Lcsr_df[Lcsr_df$obs >= Lcsr_df$lo & Lcsr_df$obs <= Lcsr_df$hi,])==0){
ggplotly(csr_plot, dynamicTicks=T) %>%
style(text = text2, traces = 4) %>%
rangeslider()
}else {
ggplotly(csr_plot, dynamicTicks=T) %>%
style(text = text2, traces = 4) %>%
style(text = text3, traces = 5) %>%
rangeslider()
}
} else if (nrow(Lcsr_df[Lcsr_df$obs < Lcsr_df$lo,])==0){
if (nrow(Lcsr_df[Lcsr_df$obs >= Lcsr_df$lo & Lcsr_df$obs <= Lcsr_df$hi,])==0){
ggplotly(csr_plot, dynamicTicks=T) %>%
style(text = text1, traces = 4) %>%
rangeslider()
} else{
ggplotly(csr_plot, dynamicTicks=T) %>%
style(text = text1, traces = 4) %>%
style(text = text3, traces = 5) %>%
rangeslider()
}
} else{
ggplotly(csr_plot, dynamicTicks=T) %>%
style(text = text1, traces = 4) %>%
style(text = text2, traces = 5) %>%
style(text = text3, traces = 6) %>%
rangeslider()
}